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Abstract. 

We study the entanglement entropy of a block of contiguous spins in excited 
states of spin chains. We consider the XY model in a transverse field and the 
XXZ Heisenberg spin-chain. For the latter, we developed a numerical application 
of algebraic Bethe Ansatz. We find two main classes of states with logarithmic and 
extensive behavior in the dimension of the block, characterized by the properties 
of excitations of the state. This behavior can be related to the locality properties 
of the Hamiltonian having a given state as ground state. We also provide several 
details of the finite size scaling. 
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1. Introduction 

The study of the entanglement in the ground-states of extended quantum systems 
became a major enterprise in recent times, mainly because of its ability in detecting 
the scaling behavior in proximity of quantum critical points (sec e.g. Rcfs. [1, 2, 3] 
as reviews). The most studied measure of entanglement is surely the entanglement 
entropy Sa, defined as follows. Let p be the density matrix of a system, which we take 
to be in the pure quantum state \^), p = |^')(^'|. Let the Hilbert space be written 
as a direct product H = Ha d) Hb- ^'s reduced density matrix is = Tr^ p. The 
entanglement entropy is the corresponding von Neumann entropy 



and analogously for Sb- When p corresponds to a pure quantum state Sa = Sb- 

The entanglement entropy is one of the best indicators of the critical properties 
of an extended quantum system when A and B are a spatial bipartition of the system. 
Well-known and fundamental examples arc critical one-dimensional systems in the 
case when A is an interval of length £ in a system of length N with periodic boundary 
conditions. In this case, the entanglement entropy follows the scaling [4, 5] 



where c is the central charge of the underlying conformal field theory and c'l a non- 
universal constant (the behavior for A'' ^ cxd is known from Refs. [6, 7]). Away from 
the critical point. 5*4 saturates to a c;onstant value [7] proportional to the logarithm of 
the correlation length [4]. This scaling allows to locate the position (where Sa diverges 
by increasing £) and a main feature (the value of the central charge c) of quantum 
critical points displaying conformal invariance. The entanglement entropy of disjoint 
intervals gives also information about other universal features of the conformal fixed 
point related to the full operator content of the theory [8] . 

Conversely, only little attention has been devoted to the entanglement properties 
of excited states (with the exception of few manuscripts [9, 10, 11, 12]), although it 
is a very natural problem. Here we consider two topical spin-chains [13] to address 
this issue. We first consider the XY model in a transverse magnetic field. We employ 
the well-known mapping of the model to free fermions to reduce the calculation of 
the entanglement entropy to that of the eigenvalues of a Toeplitz matrix on the lines 
of the ground-state case [7, 14, 15, 16, 17, 18, 19]. In the present computation, the 
properties of the excitations above the ground-state will strongly affect the form of 
the reduced density matrix and of the entanglement entropy. Then, to consider a 
truly strongly interacting quantum model, we address the same problem for the XXZ 
chain, always remaining in the realm of integrable systems. In fact, this model is 
exactly solvable by means of Bethe Ansatz [20, 21]. This provides a classification of 
all eigenstates and their energies, but no information about dynamical properties. To 
overcome this limit, we take advantage of recent progresses in the algebraic Bethe 
Ansatz [22, 23] that provides all elements of the reduced density matrix as a (huge) 
sum of determinants whose entries are func;tions of the Bethe rapidities. However, in 
this approach an inhomogeneous coupling must be considered and the homogeneous 
limit (in which we are interested) is recovered in a cumbersome manner. 

In the study of the entanglement properties of excited states, a first subtle 
point is the choice of the basis of the Hilbert space. In fact, while the ground- 
state of a local Hamiltonian is usually unique (or with a finite small degeneracy, 



Sa = -Tr pA log pA 



(1) 
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when some symmetry is not spontaneously broken), the excited states can be highly 
degenerate. Thus, any linear combination of them is still an eigenstate. In principle 
the entanglement properties can vary a lot with the basis. However, we will show that 
some of our findings are general features of all lineax combinations of the same class 
of excited states. This is not only an academic subtlety, because the exact studies 
one can perform are limited to integrable models, for which it is well-known that the 
degeneracy is large. Oppositely, any small integrability breaking term will remove 
these degenerations and one could wonder whether the specific properties found are 
only features of integrable models. 

The quantification of the entanglement in excited states can have consequences 
in the understanding of the quantum out-of-equilibrium physics and in particular of 
the dynamical problems known as quantum quenches. In fact, it has been argued that 
the post-quench state is a time-dependent superposition of eigenstates that in the 
thermodynamic limit have the same energy [24]. It is known that for a global quench, 
the entanglement entropy first increases linearly with the time and then saturates to 
a values proportional to the length of the block £ [25, 26, 27]. We will indeed find a 
full class of excited states having an extensive entanglement entropy and those could 
be the relevant ones for quench problems. Oppositely in local quantum quenches the 
asymptotic state displays a logarithmic entanglement entropy [28] and a different class 
of states should be relevant. We also mention that some of the features we find have 
similarities with what obtained in some non-equilibrium steady states [29]. 

The manuscript is organized as follows. In the next section 2 we study the XY 
model and we find two main classes of excited states, corresponding the extensive and 
logarithmic behavior of the entanglement entropy. In Sec. 3 we consider the XXZ 
model and the algebraic Bethe ansatz approach. We find that the states that have a 
logarithmic behavior in the XX limit conserve this property with the same prefactor 
of the logarithm and with a constant term slightly depending on A. Finally in Sec. 4 
we summarize our main results and discuss problems deserving further investigation. 

2. The XY model in a transverse magnetic field 

We start our analysis by considering the XY spin chain of length N with periodic 
boundary conditions, whose Hamiltonian is given by 

JV 

HxY = - 51 

1=1 

where cr" are the Pauli matrices at the site I. h is the transverse magnetic field and 7 
the anisotropy parameter. For 7=1 the Hamiltonian reduces to the Ising model, while 

for 7 = to the XX model. The diagonalization of this Hamiltonian is a standard 
textbook exercise. First a Jordan- Wigner transformation 

\m<l / \m<l / 

maps the model into a quadratic spinless free-fermion Hamiltonian (i.e. with 
anticommutation relations {cpCm} = (5/^ , {ci,Cm} = 0). After Fourier transforming 
in momentum space c/,. = cie~^^''^ /VN, the so-called Bogoliubov transformation 

bl = Uk cl + ivk c-k , bk = Uk Ck - ivk c^_k, (5) 



— — (7; (7^+1 + - 



■7 



I ""l+l 



(3) 
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makes the Hamiltonian diagonal 

N-l 

1 



H= e^l^'A-;^) > (6) 



where we considered N to be odd. We ignored a boundary term that gives a vanishing 
contribution in the thermodynamic Umit. Here we introduced the BogoUubov variables 
Uk = cos^fe/2, Vk = sin^/;/2 and angle 

^ J7sini^fc . , 27rk 

tan6'fc=- with(^fe = -— , (7) 

J cos ifk — n N 

giving single-particle eigenvalues 

\J{h- J cos LpkY + sin^ ■ (8) 

From this dispersion relation, it is evident that the model is critical (gapless) for 7 = 
and \h\ < \ J\ (XX universality class) and for h = ± J and any 7^0 (Ising universality 
class) . 

The exact diagonalization of the model gives not only the ground-state properties 
but a complete classification of all the eigenstates and in particular their energy. In 

the basis of free fcrmions, the excited states arc classified according to the occupation 
numbers of the single-particle basis (that is the basis of Slater determinants). A 
generic eigenstate can be written as 

= n ^llO) ' ^itli energy Ee^ = U Sk - ek \ , (9) 

where is the set of occupied momenta. To give a simple pictorial representation 
of these states, we indicate with up-arrows the occupied single-particle levels 
(excited quasiparticles) and with down-arrows the empty ones, with the first arrow 
corresponding to momentum ipk = — tt. When a set of n consecutive momenta are 
occupied (empty), we simply replace the up (down) string with |" (i")- For example, 
the ground state is | J, . . . J.) = | J,^). Counting all the possible arrow orientations, 
it is obvious that this graphical representation generates all the 2^ eigenstates of the 
chain. Notice that these arrows have nothing to do with the state of the spin in real 
space (the real space configuration is an highly entangled superposition). 

When calculating the entanglement entropy, three different length scales enter 
in the computation: the size of the chain N, the length of the block £ and the 
number of excited quasiparticles that is encoded in the size \Ex\ of the set E^. 
General results can be obtained in the thermodynamic limit AT — > 00 and when 
£ » 1 (in finite size, this limit describes the regime N ^ £ ^ 1). In this 
limit, it is obvious that if only a small number of quasiparticle levels are populated 
(i.e. \Ex\ <C N), the corrections to the ground-state correlation matrix can be 
generally treated perturbatively and in a first approximation the excited quasiparticles 
contribute independently to the entanglement, giving rise to a negligible contribution 
in 1/N. Thus, in the thermodynamic limit, all entanglement properties of these states 
arc equivalent to those of the ground- state, but this does not prevent from interesting 
and maybe calculable finite-size behavior. We have been informed of some unpublished 
work by M. Ibanez and G. Sierra [30] studying the entanglement entropy of these low- 
lying excited states obtaining a finite-size scaling different from Eq. (2). Here instead 
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we are interested in those states that are macroscopically different from the ground- 
state and that will have an entanglement entropy that in the thermodjmamic Umit 
could differ strongly from Eq. (2). 

In order to work directly in the thermodynamic limit, we need a proper description 
of excited states. This is rather straightforward. In fact, when N ^ oo the possible 
values of k are all the integer numbers and the reduced momentum ipk becomes a 
continuous variable (/? living in the interval <f G]—Tr, 7r[. We are here interested in the 
case with jE'^.j ^ N (that can be seen as an "highly excited state", even if it is not the 
energy that matters). Thus in all the formulas involving sums over populated energy 
levels, we substitute sums with integrals by using as distribution a proper defined 
regularized characteristic function of the set that we will indicate as m(ip). The 
function (1 + m{ip))/2 represents the average occupation of levels in an infinitesimal 
shell around the momentum ipi- = 2-Kk/N. Let us give several examples to make this 
limiting procedure clear (a < 1): 

I i"") m(^) = -1 , 



1 , if < 9? < Tra , 

— 1 otherwise , 



1—1 otherwise, 

|{Ti}^/') 

otherwise. 

We only wrote down for simplicity states with a step-wise characteristic function, but 
with little fantasy it is easy to imagine states with a smooth one |. 




2.1. The reduced density matrix and the entanglement entropy 

It has been shown [7, 14] that, despite the non-local character of the Jordan- Wigner 

transformation, the spectrum of the reduced density matrix pA of a single interval 
A = [0,£] is the same in the spin variables cr; and in the free-fermion ones c;. This 
property makes the XY model the ideal testing-ground to understand the behavior of 
the single-block entanglement for the excited states. The eigenvalues of the reduced 
density matrix pt of a block of £ adjacent spins for a Slater determinant are related to 
the eigenvalues Vi of the correlation matrix restricted to the subsystem [7, 14] . This is 
easier to sec by introducing the Majorana operators Af = cj + and Af = i{ci — cj) 
[7]. The eigenvalues of pi can be labelled with the configurations of £ classical spin- 
variables denoted as Tj = ±1 and it holds A{^} = 11^=1(1 + 'O'^j)/^; with ii'j the 

f If we would be pedantic in defining tfiis limit, we can think to (f + m{ip))/2 as the convolution of 
the characteristic function of Ex with a Gaussian of zero mean and standard deviation that must be 
put to zero at the end of any computation. Since in the sum in Eq. (12) there is almost everywhere 
(everywhere in non-critical regions) a regular function of <p, the regularization in the definition of 
m(ip) is perfectly well-defined. 



Entanglement entropy of excited states 
eigenvalues of the block Toeplitz matrix 

/To ••• r^_i \ 

n = 



A- 



)( (11) 



\ ri_^ • • • To / 

The two-by-two matrices F; arc easily computed observing that the generic eigenstate 
in the Slater-determinant basis (9) is the vacuum of the fermionic operators 

k e E^, 

^ h\ otherwise. 

After simple algebra one obtains 

p(GS) 2i ^ / sm{lipk) - cos{lipk - Ok) \ 
' ~ ' ^ N /^\cos{lvk + ek) sin(/^fe) J' ^'^> 

where 9k is the Bogolioubov angle of the transformation that diagonalizes the 
Hamiltonian in Eq. (7) and r^^'^'^-' the corresponding matrix in the ground-state [7, 14] . 

As explained in the previous subsection, when \Ex\ ~ A'', we can substitute in 
equation (12) the sum with an integral 



1 r , 1 + m{<p) 



N ^ 27r 



d</3 ^ iPk^^, (13) 



where (1 + ■m{(f)) /2 is the regularized characteristic function of the set E^ introduced 
above. Substituting in Eq. (12) this regularization we have 



^J_r d^e-*'^r(f-) , with 



(14) 



p(£V) ^ 1 ■m{-ip) - ■m{ip) -i[m{ip) + m{-ip)]e''^ \ 

The entanglement entropy can be expressed as a complex integration over a 
contour C that encircles the segment [—1,1] at the infinitesimal distance r/ as in Ref. 
[17] 



where 



\ f d 
Se^ lim -— * dAe(l + 277,A)-rT-logdet|Al-n| 

7]^o+ Am Jc dA 



, , X + y , x+y x—y, x—y 

e(x, y) = log log 

\ yyj 2 2 2 2 



(16) 



A similar expression is easily written for all Renyi entropies for general n. Applying 

the Szcgo lemma (see e.g. Ref. [31]) to the determinant of the block Toeplitz matrix 
Al — n, we obtain the leading order in £ of the entanglement entropy 

St=^f H{m{v)) + 0(log I) , (17) 

with H{x) = e(l,,x-). 

This first result is very suggestive: the entanglement entropy of a class of excited 
states in the XY model is extensive, in contrast with the logarithmic behavior of the 
ground state. However, every time that m((p)'^ ^ 1 only in a region of vanishing 
measure of the domain (as in the ground state) this leading term vanishes, and one 
should go beyond the Szego lemma to derive the first non-vanishing order of the 
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entanglement entropy. It is important to stress that for this type of "highly excited 
states" the leading order of the entanglement entropy is not sensitive of the criticality 
of the ground-state. This does not come unexpected, because we are exploring a region 
of energy that lies extensively above the ground-state. 

To describe the (sublcading) logarithmic terms in the determinant of a Toeplitz 
matrix, we should use the so-called Fisher-Hartwig conjecture [33]. If m(<p)^ = 1 
almost everywhere, m{ip) can be re-written in the following form, that is particularly 
useful to apply Fisher-Hartwig e] — n, 7r[) 

n 

m((^) = e'*'^SmW'Qgiarg(v=-v:.^)^ j-j^g) 
J=l 

where 2[n/2]§ is the number of the discontinuities of m(ip) and ipj are the discontinuity 
points (the term 2\n/2\ takes into account an eventual discontinuity in tt that is not 
counted by considering the open interval (p g] — n , Tr[) . We prove analytically in the 
next subsection that Sg oc log £ in the XX chain (7 = 0) and then we show that this 
is not a peculiar feature of the isotropic model. 



2.2. XX chain 

In the XX spin chain the Bogolioubov angle reduces to e*^*^ = sign(Jcos(^fe — h) 
and the Fisher-Hartwig conjecture is sufficient to prove the following result: the 
entanglement entropy of the excited states described by the multi-step function (18) 
grows logarithmically with the width of the block. The coefficient in front of the 
logarithm is 1/6 times the number of discontinuities in the non-critical region (|/i| > 1) 
and it must be corrected in the critical region (|/i| < 1) to take into account the modes 
with zero energy. For \h\ < 1 the modes with zero energy at ±(pF {Vf = arccos \h/J\) 
define the function 

= < . . . (19) 

I —m{—(p), otherwise, 

that substitutes m{(f) when counting discontinuities. The importance of the number 

of discontinuities was firstly stressed in [32] in a different context. In Fig. 1 a direct 
computation shows the importance of the position of the modes with zero energy. 



2.2.1. Fisher- HaHwig proof of the log-behavior in the XX chain. The proof of the 
relation between the entanglement entropy and the discontinuities of m{ip) when Eq. 
(18) holds (i.e. when m{(f) = ±1) in an XX chain is a slight modification of the proof 
given by Jin and Korepin in Ref. [17] for a critical XX ground-state. For 7 = 0, the 
matrix (15) can be written in terms of the Pauli matrix ay as 

T{ip) = ±aym{^cryip) , (20) 

with the upper (lower) sign if the momentum if is below (above) the Fermi level of 
the Jordan- Wigner fermions. As a consequence the block Toeplitz matrix (11) can be 
reduced to a standard Toeplitz matrix with symbol 

^^^^^fl, (e^>OAm(-<^) = l) v(e^<OAm(<p) = -l), ^^^^ 

1—1, otherwise. 



§ Here and below, \x] stands for the closest integer larger than x and [x] for the closest integer 
smaller than x. 
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Figure 1. The entanglement entropy as a function of tiie block length for the 
excited state with characteristic function m{tp) = sign(((/3 — -g)!^ — of two 
critical XX chains. The different behavior is caused by the position of the zero 
modes {ipp = 7r/6 and tt/2 with two and four discontinuities respectively) and 
results in a = 1 or 2. The straight lines are the analytic prediction for large £ 
given by Eqs. (25) and (26). 



with e^p — J cos (fi — h. The reduced correlations matrix Al — 11 is generated by the 
symbol 

n 

where the (^j's are the momenta corresponding to the n discontinuities of 7(</3). The 
ground state has two symmetric discontinuities at ±</?i?. The symbol admits the 
canonical Fisher-Hartwig factorization [33] 

n 

i(^) = (A+ir(A-l)''ni.M. 



with 



(fij < (fi < ipj + 27r 

rA + l 



2ni 



X-1 



TT < arg 



A - 1 



(22) 
(23) 



and the two exponents are 



n 



Defining kp = ^fj/'^^ ^^e Fisher-Hartwig conjecture (that for this 

case with |A| > 1, i.e. |Re(/3j)| < 1/2, has been proved by Basor [33]) reads 



det |A1 



■n 

ni-n[ 



(2 - 2cos((^i - ipj)) 



(24) 
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xG(l + /3(A))"G(1 -/3(A))"{(A + '"^"| V"'^^^)' 

where = /?? and G{x) is the Barnes G-function 

G(l + /3)"G(1 - = e-(i+^«)"''' n - ^) '"e"^'/'' • 

In order to find the entanglement entropy we have to evaluate g^log-D£(A), where 
De{X) = det |A1 — n| (cf. Eq. (16)). The derivative can be easily computed and it 
consists (in principle) of three terms giving in Eq. (16) 

Se^aoi+^logi+a^^^}, (25) 

with: 

- the linear term aoi is the same as in the ground state [17] (except from the 
definition of kp), and it is known to vanish ao = (as actually we already proved); 

- the logarithmic term a/3 log £ is the ground state contribution multiplied by 
a = n/2 (a will be interpreted as an effective central charge, that is why we 
multiplied by 1/3); 

- the additive constant a^^^^ is slightly more complicated but it has essentially the 
same structure of the ground-state and it is 

n -^(-1)^"*, \ . 2^'Pi-'Pj\^ /o«^ 
= 2 ~ 2^ 6 [sm ( ^)J , (26) 

i<j 

and ao is the additive constant for the entanglement entropy of the critical 
XX chain without magnetic field ao ~ 0.726 . . . (see Ref. [17] for the analytic 
expression) . Notice that it depends not only on the number of discontinuities but 

also on their location. 

See Fig. 1 for a comparison of this analytic asymptotic result with the direct 
computation for finite i. 

In Ref. [19] Igloi and Juhasz showed that the ground-state entropy of the XY 
model with h = can be related to the sum of two Ising models (i.e. 7 = 1) with 
fields h depending on 7. When specialized to the XX model, the two Ising chains are 
both critical and one has 

Sxx {2£, 2N) = 2Sisinsie, N) . (27) 

The proof in Ref. [19] can be generalized to some excited states by properly rescaling 
all length scales. Thus the knowledge of the result for the XX model, automatically 
gives the value for the critical Ising chain. However, we do not show the details of 
this proof here, because in the following we will provide the asymptotic result for any 
logarithmic state of the XY chain. 

At this point it is natural to wonder whether these eigenstates having an 
entanglement entropy growing logarithmically with £ are the ground-states of some 
conformal Hamiltonians. In the case of the XX model, since Hxx{h) with different 
magnetic fields commute among each other, the ground-state at given h is an excited 
state of a chain at different h. Thus, for all these states it is obvious that they should 
display an entanglement entropy scaling like Eq. (2) with c = a = 1, i.e. they have 
two discontinuities in m{ip). As we will see, this is true in general and in the next 



Entanglement entropy of excited states 



10 




Figure 2. The entanglement entropy as a function of tlie block length for two 
excited states of the infinite critical Ising chain with 4 discontinuities (Left) at 
momenta {0, 0.5, 0.8, 1.4} and {—0.5, 0.5, 0.8, 1.4}. The different slopes are caused 
by the zero mode. Right: Two excited states of a non-critical XY chain in finite 
size. 



subsection we show that a commuting set of local operators of the XY chain can be 
used to prove that all these logarithmic excited states are ground-states of properly 
defined local conformal Hamiltonians. Eq. (25) can be exploited to deduce the central 
charge of this local Hamiltonian c — a — n/2. 

2.3. Logarithmic behavior and effective Hamiltonians 

It is straightforward from Eq. (12) to calculate the spectrum of the reduced density 
matrix and the entanglement entropy for any eigenstate at any value of 7 and h. 
We calculated the entanglement entropy numerically for several different cases and we 
always find a logarithmic behavior with ^ every time m{iff' = 1 almost everywhere (see 
e.g. Fig. 2). To get a proof similar to the one of the previous section for the general 
XY model, one should generalize the methods in Ref. [18] mapping the computation 
to a Riemann-Hilbert problem. This way of proceeding is very complicated and 
we take here a different route based on the considerations we reported at the end 
of last subsection. In fact, this general logarithmic behavior of the entanglement 
entropy suggests that this type of excited states can be the ground states of critical 
Hamiltonians. We explicitly build these critical, translational invariant, and local 
Hamiltonians, proving the logarithmic behavior, with the correct prefactor. 

The excited state \Ex) in Eq. (9) is the ground state of all free-fermionic 
Hamiltonians of the form 

H = Y1 e{ipk)bibk, with e{ipk) <0<^keE^, (28) 
fc 

for any choice of the function s{ipk). In particular we could choose s{ipk) — 
—f{ipk)m{Lpk), with f{x) an arbitrary positive function. The choice of s{ipk) 
determines the locality properties of H: most of the choices of e{(pi~) would produce a 
non local H (while by construction H is always hermitian and translational invariant 
because it is built by Fourier transform). 

To understand the locality of this effective Hamiltonian it is useful to introduce 
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the operators (^^'^ are the Majorana operators introduced above from Ref. [7]) 

G(r) = i ^ AfAl^ , and F^(y^ (r) = * ^ Af^^A^ll^ . 
I I 

In fact, by separating e{ipk) in its even and odd part (e(<Pfe) = ed'-Pk) + Soifk)), we 
can rewrite the effective Hamiltonian as the sum H = + Ho where 

iV-l 

r h— ^ — ^ r 

N-l 

^° = ^E[^ E eMe-"''''-]{F:: + Fy)^Y.9o{r){F: + Fy), (29) 

r 1 i-w r 

n— 2 

where we defined the complex coupUngs ge{T) and Qoix). 

The locality of H is related to the long distance behavior of these complex 
couplings ge/oii")- From a standard theorem in complex analysis, wc know that 
Qe/oif) decay faster than any power (and so results in local couplings) if their Fourier 
transforms in the above equations are C°° (i.e. with all derivatives being continuous 
functions; often wc will refer to these functions simply as regular). When Eq. (18) 
holds, that is m{ip) = ±1 has a finite number of discontinuities, and for a non-critical 
system (i.e. when e"'^*^ is regular), the arbitrariness in the choice of e allows us to 
take it among the functions. This conclude the proof for non-critical systems. 

For the critical case, a slight modification is enough to give the correct 
Hamiltonian. In the XX spin chain e~'^ = sign(Jcos(/? — h) so that we can make 
the two above functions regular simply defining the characteristic function m(y) 



Th{(p) 



—m{—ip) otherwise , 



as wc have already done in Eq. (19). The critical XY {\h\ = 1) is more involved 
because e^*^ can be made regular only after imposing anti-periodic conditions to the 
mode of zero energy. It is then convenient to extend the definition of e to the interval 
[0,47r] 

e(4^)(</') = S ^ r-to A ^ (^^) 

I — e(47r — ip) tp e [2tt, i-K] . 

e(47r) can be chosen C°° because it has at most 2n + 2 zeros, where n is the number 
of the discontinuities corresponding to the excited state. The constructed function 
restricted to [0,27r] has the correct regularity properties. Regardless of the presence 
of a discontinuity in (/? = the dispersion law must vanish in = ( see Eq. (30)), 
thus the number of chiral modes is the number of discontinuities, plus 1 if there is not 
a discontinuity in ip = 0. This ends the construction of the local Hamiltonian for all 
the XY models. 

And this is not yet the end of the story. We can in fact use the arbitrariness we 
have in the choice of ik to fix it in such a way that it crosses the zero-energy line with 
a non-vanishing slope. The low-energy properties of the resulting Hamiltonian can 
be then studied by linearizing the dispersion relation close to the zeros in a canonical 
manner. Each zero gives a chiral mode with central charge 1/2 and so the total central 
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Figure 3. Two 5-folded wrapped chains of 60 spins. The thick green line 
represents the subsystem (6 spins on the left and 18 spins on the right) while 
the red links give weight to the interaction between the subsystem and the rest 
of the chain. If the "area law" holds the entanglement entropy is proportional to 
the number of the links. 



charge will be n/2, with n the number of zeros, i.e. the number of discontinuities of 
m{ip) for non-critical systems, or the proper variation for critical ones (when the 
zero mode gives one additional contribution). This agrees with all the specific cases 
in the previous section. In particular if m(ip) is discontinuous in (/3 = 0, the zero 
mode contributes only once. In Fig. 2 we report some specific examples stressing the 
importance of the critical modes and of the location of discontinuities. 

2.4- Finite size scaling 

When the width of the block £ is comparable with the length of the chain N, the 
characterization of the entanglement becomes tricky. When an excited state \Ex) 
can be associated to the ground-state of a local Hamiltonian H with central charge 
a — n/2, i.e. when the entropy grows logarithmically with £ with a pref actor given 
by a, the constructive proof of previous subsection in the thermodynamic limit is still 
valid. Thus, in this case, the entanglement entropy has the finite size scaling given by 
Eq. (2) with c replaced by a. This is shown in the right panel of Fig. 2. 

A more intriguing problem is to understand the finite size scaling of excited states 
that have an extensive entanglement entropy in the thermodynamic limit. The result 
for N —I- oo only predicts the derivative of the entropy for small subsystems. Increasing 
£ peculiar finite size behaviors must emerge, because the chain is finite and the entropy 
must be symmetric around £ — N/2. 

Up to now we studied in detail excited states with a regularized characteristic 
function of the type (18), that is \Ex) = \ 0^=1 T"n™^)i where rij and mj are all 0{N) 
and d is a finite number. States with m(ip)'^ ^ 1 (that have extensive entanglement 
entropy) do not fall in this category as evident in the definition (10). They can be 
realized by joining in a regular fashion small blocks k made by a given sequence of 
populated or empty energy levels (e.g. k = {ti} or k = {t^i} etc.). Thus to study 
the finite size scaling of "extensive" states, we concentrate on those of the form 

d 

\E,) = \l[n"^R'^^), (31) 
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where R is the set obtained interchanging f with J.. The entanglement entropy of 
this type of states in the thermodynamic limit has an extensive behavior because k 
averages to give m{ip) = {u — d)/{u+d), where u {d) is the number of up (down) arrows 
in K, while R gives m{if) = [d — u) / {u + d): the regularized characteristic function is a 
multi-step function but with modulus different from 1 and Eq. (17) gives the leading 
term of the entanglement entropy. 

In order to have a quantitative prediction for the finite size scaling, we follow the 
ideas in the previous subsection by looking at the effective Hamiltonian obtained by 
the construction in Eq. (28). The resulting couplings in Eq. (29) could never give 
a finite-range Hamiltonian because the entanglement entropy is not logarithmic. We 
can make a local choice of the sign that makes i a regular function (that we call e) 
giving the coupling || 



E ^"""'^i.^k) (32) 

1 r 1 1 



N 



and the interaction is not local anymore. The 0{1/N) term comes from the series 
expansion of e. The first factor in equation (32) is periodic of period N /\k\ while the 
second one is a modulation. The coupling decays faster than any power for r < A^/2|k|, 
but it explodes (i.e. it grows faster than a power) up to N/\k\ when g{r) becomes 
again of order 1. The behavior for large distances is determined only by the first region 

The interaction is localized within a distance iN/\K\. so the Hamiltonian can be 
interpreted as a local one in a |/t|-folded wrapped ID chain. If we assume the "area law" 
to be valid for the wrapped chain (i.e. that only a shell of mutually interacting spins 
contributes to the entanglement [34]), we can predict the behavior of the entanglement 
entropy: each spin strongly interacts with the neighborhood spins and with the \k\ 
spins of the other wrappings (see Figure 3). Thus the entanglement entropy is a 
piece-wise function of (, that changes slope at jN/\K,\. We can find excited states with 
analogous properties considering any finite partition of unity of the circle ] — tt, tt] , 
with the property that all functions of the set are regular and approach step functions 
in the limit of large N. We associate a small block k*^*^ to each function of the set and 
we write the coupling as a sum of terms of the form (32) 

n 

that we obtain identifying the regularized e with the given function of the partition. 

In the scaling limit the characteristic function is m ^ Wt^i'^l)- Each g^, has the 
behavior previously described, so the entanglement entropy is a piece- wise function 

II This coupling is a slight modification of the ones in Eq. (29). It has the advantage to make all 
the formulas simpler, but it applies only to non-critical systems. The generalization to critical ones 
is straightforward, but long and we do not report it here for clarity. However all results (except for 
the ground state) are independent of this choice, as in the previous section. 
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Figure 4. Examples of 3- and 4-folded states. Left: The 3-folded excited state 
I J.'** {TJ.'^}^^ J.*^) for tlie non critical chain {h = 0.5,7 = 0.5). Entropy grows 
linearly up to Af/3 and then saturates. The dashed line has the slope given by Eq. 
(17) with the regularized step-function m{(p). Right: The 4-folded excited state 
I J.^° {Ti^}"*^^ J.^"*^)- For ^ < ■'V/2, the entropy always grows Uncarly, but with a 
change of slope close to £ ^ 7V/4. 



of £ changing slope in where \k\ is the least common multiple of the {I'tl*-*'}- 

Two examples of 3-folded and 4-folded states are reported in Fig. 4. 

To give the details of a specific example, we report the 3-folded case k^^^ = {t^i} 
and K^^-' = {1} with (p € Ii ^ cos(p > 1/2 and ip £ Iq <^ cos<p < 1/2, in other 
words the set is made of the quasiparticles with momenta (27r(3fc + q))/N with 
|A:| < iV/12 and q e {0,1} 



1^-)= n 



The excited state \Ex) is the ground state of the Hamiltonian 

JV-l 

H= ± [(^-cos^.)(-l)rfU(-l)LfJ+(-l) 



and if N is divisible by 3 the coupling is different from only in 9 points 

r = 0, 
r = ±l. 



9{r) 



± 



1 



12 



3 



{-1,0,1} 



The effective Hamiltonian H is local on the 3-folded wrapped chain. The entropy grows 
linearly with the width of the block up to £ = iV/3, after that the interaction surface 
does not further increase and the entanglement entropy does not depend anymore on 
the width of the block, see Figure 4 (left). Notice on the same figure (right), the 
change of slope in the 4-folded case. 
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Figure 5. Rescaled half-chain entanglement entropy t = {N — l)/2 for a critieal 
XX in zero magnetic field, a critical Ising, a non critical XY spin chain, and 
an Ising in a very large magnetic field. All plots are for N = 15. Each point 
corresponds to an excited state with energy (in unit of J) on the real axis. The 
red curves are the "2-folded" estimations of the envelope. 



2.5. Some general properties 

To have a general picture of the scahng of the entanglement for all excited states and 
not only in the particular classes considered so far, we study here the entanglement 
entropy in a small enough chain to be able to calculate it for all the 2^ states. We 
mainly concentrate on blocks with maximal entropy, i.e. with length equal to half- 
chain (actually {N — l)/2 spin, because we use N odd). Drawing general conclusions 
in an analytic manner for finite systems is not easy, so we mainly analyze numerical 
results. The plots in Fig. 5 suggest that some regularities are general features of 
excited states and not only of the classes we can compute analytically. In these plots 
(and in all those relative to this section) we always consider the rescaled entropy 

rescaled entropy = , with Sf^ = - log( — sin — ^) , (33) 

SO that, for states with a critical- like behavior (for large enough £ and N) we have a 
direct estimation of the effective central charge. We found particular instructive to 
plot the (rescaled) entanglement entropy as function of the energy of the eigenstates. 
In Fig. 5, we considered chains of 15 spins and we plot the rescaled 6*7 for all the 2^^ 
eigenstates. Similar plots can be done as function of total momentum instead of the 
energy. 
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rescaled entropy 



Figure 6. Histograms for tiie number of thie states witli a given entanglement 
entropy for a non-critical XY chain of 23 spins, after cutting the Hilbert space in 
an energy shell. Main plot; rescaled 5ii . Inset; rescaled Sg. The band-structure 
is evident only for £ = 11. 

A first feature that is particularly evident from the plots is the band-like structure 
of the entanglement entropy (notice that this is independent of the use of the energy on 
the horizontal axis, any other conserved quantity would result in qualitative similar 
plots). This means that the entanglement entropy of excited states distributes at 
roughly integer (or half-integer for critical XY at h — 1) multiples of Sf^. For states 
with a small number of discontinuities (compared to N) , this phenomenon is clearly 
due to the quantization of the prefactor of the logarithm. However, in general this band 
structure cannot be so easily explained: the excited states with a logarithmic behavior 
are expected to be negligible in number compared to all the others. Increasing the 
number of discontinuities at fixed N, the crossover to extensive behavior takes place 
and eventually it deteriorates the bands. This last phenomenon is not evident in 
Fig. 5 because the band structure persists up to the maximum allowed number of 
discontinuities. The simplest explanation is that also extensive states should roughly 
be quantized but within a scale different from Sf^, that in particular does not grow 
with N. To check this, we should increase N, but in doing so, the dimension of 
the Hilbert space grows exponentially and it becomes soon prohibitive to plot (and 
understand) so many points in an readable graph. For this reason we considered a non- 
critical chain of 23 spins, and to reduce the number of states, we limited to states with 
energy in the interval 4.600 < E ~ Eq < 4.694. In Fig. 6 we report the distribution of 
the points. For £ = 11, the band structure is evident and the points distribute in an 
almost Gaussian fashion around some discrete values of the entanglement entropy, but 
the distance between them becomes smaller Sf^, confirming that the origin of this 
phenomenon in the upper part of the band has nothing to do with logarithmic states. 
For i = 6 (inset of Fig. 6) the band structure disappears completely, confirming that 
most of the states are extensive. We checked that still increasing N, this scenario is 
consistent. 

Another very interesting feature is that in all the plots, the entanglement entropy 
has a maximum value that seems to be a regular function of the energy (that is the 
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energy energy 

Figure 7. Rescaled entanglement entropy for small blocks. Left: £ = 4 in a 
non-critical XY-chain of 15 spins; The continuous curve is Eq. (34) giving a good 
estimation of the envelope. Right: £ = 5 in a non-critical Ising chain of 15 spins; 
The "3-folded" envelop (in red) of the envelope is in good agreement with the 
data. For high energies, when the "3-folded" approximation is not defined, Eq. 
(34) (in green) works well. 

final reason why we made this kind of plots) . We argue here that these envelopes have 
a characteristic dependence on the energy that in the scaling limit is determined 
by excited states with extensive behavior. We already derived the entanglement 
entropy for the excited states that are equivalent to the ground state of n-folded 
wrapped Hamiltonians. Eq. (17) characterizes the scaling regime, e.g. for the 2- 
folded case the entanglement entropy increases linearly up to N/2, while in the 3- 
folded it increases up to N/3 and then saturates. We have then for blocks of length 
£/N > iJ(l/3)/(27?(0)) = 0.459 . . . that the 2-folded case is more entangled than the 
3-folded one. This suggests that the 2-folded states can explain the envelopes in Fig. 
(5) for £ = {N — l)/2 (Notice that the maximal entangled state, regardless of the 
energy, is always a 2-folded one). If this is true, the envelope is easily obtained: the 
problem is analogous to find the dependence of the particles density on the Fermi 
energy in a free Fermi system at zero temperature. Indeed using Eq. (17) and the 
asymptotic expression for the energy, the "2-folded approximation" of the envelope 
satisfies the parametric equations (valid for E < 0) 

( 5*^^^ log 2 r 

[|--^£d^^^(e-M). 

In Fig. 5 this analytical result is compared with the numerical data for a critical XX, 
a critical Ising and two non critical XY spin chains: the approximated envelope is in 
good agreement with the numerical data also for small chains. We also notice that 

dSi^ ^ log2 

dE - ~ V dE Jg.s. ' 
where A is the gap in the dispersion law: if the system is critical then the "2-folded" 
approximation of the envelope has infinite derivative in E = Eq g (cf. Fig. 5). 

In the opposite limit of small £, the band structure is practically lost (see left 
panel of Fig. 7) and for most of the states Eq. (17) gives a good estimate of St so 
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that we expect that the envelope can be determined maximizing the expression (17) 
at fixed energy. The maximization gives the thermal-Uke parametric equations 

gMAX ^ r-K 

dip il(tanh(/3e)) , 

(34) 

dip etanh(/3£) , 



t 2tt 

E 

N ' 47r 



-ir, 



and the loss of the band structure can be seen as a consequence of a "pure" extensive 
behavior of the entropy. Eq. (34), in the scahng hmit, is always an upper bound for 
the entanglement entropy because entropy is a concave function of i. In Fig. 7 (left) 
we compare this analytical curve with the data for AT = 15 and ^ = 4 in a non-critical 
XY-chain. 

Considering blocks of intermediate lengths the parametric equations (34) define 
a too high bound (see right of Fig. 7). At the same time the band structure starts 
emerging. We can improve our estimation considering a generalization of the "2-folded 
approximation" of the envelope; the "n-folded approximation" (that makes sense only 
for € < N/n). The maximal entanglement entropy in the folded family of excited 
states is 



gMAX H{1 - I] 



f 

J — TT 



dip e{n - s) , 



N 2nn 
N 2mr 



l^£d, -5). 



In Fig. 7 (right) we report ^5 for a non-critical Ising chain of 15 spins (so the maximum 
allowed n is 3). It is evident that up to the point where it exists the 3- folded curve is 
a good approximation of the actual envelope, while for larger values Eq. (34) works 
well. 

All the plots in this subsection are relative to the Slater-determinant basis. 
We have checked that considering linear combinations of eigenstates with the same 
energies, these envelopes remain unchanged, while the band-structure disappears (as 
maybe could have been expected). 

Lack of space prevents us to show many other similar plots about the distribution 
in the energy of excited states for the entanglement entropy. The main features 
about appearance and disappearance of the band-structure and the envelopes (that 
we showed here with few examples) are always true. 



3. The XXZ model and the algebraic Bethe Ansatz approach to reduced 
density matrices 

We consider the anisotropic spin- 1/2 XXZ model in the presence of a magnetic field 
in the z direction, with Hamiltonian 

N ^ 

fez = {«+! + «+i + - 1) - 2<}' (36) 

m=l 

and periodic boundary c;onditions. The model is solvable by means of the Bethe 
Ansatz for any real value of the anisotropy parameter A [21, 35], but we will consider 
here only the antiferromagnetic critical regime < A < 1 (the case A = is the XX 
model of the previous section). We will use the quantum inverse scattering solution for 
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this problem found by Kitanine, Maillet and Terras [22, 23]. Some further advances 

for the algebraic Bethc Ansatz approach to the XXZ model can be found in [36]. It 
is worth mentioning that recently also the full solution for open boundary conditions 
has been found [37]. 

The approach we follow in this paper is inspired to the ABACUS method 
pioneered by J.-S. Caux and collaborators to calculate exact dynamical correlation 
functions in finite systems [38, 39, 40, 41, 42, 43]. In fact, instead of searching 
for exact relations valid in the thermodynamic limit (as e.g. made in few specific 
cases for the ground-state properties [44, 45, 46, 47]), we will work with finite chains, 
solve numerically the Bethe equations for a given eigenstate and plug the solutions 
in the determinant form found for the elements of the reduced density matrix. The 
computation of the final result will require many computational resources and we 
will discuss in the specific cases when this approach is more convenient than exact 
diagonalization. We mention that for A = 1/2 and A'^ odd, thanks to very peculiar 
combinatorial properties [48], some exact results are known also for finite chains 
[49, 50]. 

The content of next subsections is highly technical. We first review the main 
results of Ref. [22, 23] (to make this paper self-consistent and to fix the notations) 
and then we explain the technical tricks to adapt these fundamental results to the 
computation of the reduced density matrix. We remand the reader interested only in 
the results to the final subsection 3.5. 



3.1. The algebraic Bethe Ansatz 

In the algebraic Bethe Ansatz approach (see the book [35] for an introduction to the 
subject), the dynamics of the model is encoded in the so called R matrix 

/I \ 

b{X,^i) c(A,/i) 

c(A./i) 6(A,/7,) 

\ 1 J 



R{X,li) 



(37) 



where 



6(A,/x) = 



sinh(A — /i) 



sinh(A — fi + T])' 
Here the parameter 77 is related to A by the relation 

A^^-{e^ + e-^). 
Now we introduce the monodromy matrix 

T(A) = Ron{X - ^jv) ■ . ■ -Ro2(A - 6)^oi(A - ^1) 



c{X,fi) = 



sinhr/ 



sinh(A — n + rj)' 



(38) 



A{X) B{X) 
C{X) D{X) 

where are arbitrary parameters sitting on each site of the spin chain. The role of 
the inhomogeneities will become clear in the following. We introduce the transfer 
matrix as trace of the monodromy matrix T(A) = TrT(A) that satisfies 



[limr(A,0. lim T(A,0] =0, 



with a= {a,. . . ,a) , 



(39) 



where we denoted with ^ the vector with components In the approach of Ref. [22] 
keeping the different helps in deriving general results. Only at the end, to recover 
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the results for the XXZ model, we will take the so called homogeneous limit ^ — > a. 
Every eigenstate of the Hamiltonian (36) can be written as 

M M 

I {Ai} > = n B{X,)\ ), ( {AO I = ( I n C(Afe), (40) 
fe=i fe=i 
where we denoted with 1 ) the reference state with all spins up 

N 

|0) = (g)|+)fc. (41) 

k = l 

The parameter M is such that M < N/2 and M = N/2 for the ground state in 
zero magnetic field, while the parameters {Ai, . . . , Am} are called rapidities. We also 
introduce d{X) 

N 

d{X) = Y[ b{X, Ci) , for which d{^i) = Vi . (42) 
1=1 

Not all states of the form Eq. (40) are eigenstates of the Heisenberg Hamiltonian: the 
rapidities Aj must satisfy a set of non-linear equation known as Bethe equations that 
for the Heisenberg chain can be written as 

1 |^6(A^^ _ 
^(Aj) t=i KAfc,Aj) 

We need the commutation relations 

[B{X),B{n)] = [C(A),C(m)] = 0, for all A,m, (44) 
to derive the action of the operators A,B ,C ,D on a.n arbitrary state | {A^} ) [23] 

M 

M M+1 n sinh(Afe - Xa' + rj) m+i 

(0|nC(Afe)^(AM+i) = E «(^«')'^^^^ (0| n C(Afe); (45) 

'==1 n sinh(Afc-A„0 

Je = l 

M 

M M+i n sinh(Aa -Xk+rj) m+i 
(0|nC(Afe)^(AM+i)= E^(^«)^^b ^'^1 n (46) 



^==1 «=i n sinh(Aa-Afe) 



fc=i 

k = -L 
k^a 



M 



M M+1 n sinh(Aa -Xk + vi) 

(0|nC(Afe)S(AM+i)= E'^(^»)^1^T^ 

^==1 «=i n sinh(A„-Afc) 



k = l 
k^a 

M + 1 

n sinh(Aj - Xa' + r?) 

M+1 . /Jl M + 1 

i: .„mam°'-v , — -(01 n^^')^ m 

M+1 ^ sinh(Aj-A.O 

i = i 



We can fix a(A) = 1 for all A. 
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A fundamental ingredient is the formula for the scalar product between two 

arbitrary states. Given a set {Ai,...,Am} that is solution to the Bethe equations 
(43) and another set of arbitrary numbers {/xi, . . . ,/xm}, the scalar product of states 
of the form (40) is given by the so called Slavnov formula [51] 

detiJ({AJ,{/i,}) 



(0, n-te) u^iKm - n...,:-;:;!^^^;.,-.. 

^~ j>k a<0 



(48) 



0' 

j>k a</3 

where we defined 

Hab = . l^!^!^^^'' X I :77^ IT sinh(A„j -fXh + rj)- TT sinh(A„j - /Xft - ??) ) . (49) 
smh(A„ - fib) \ d{nb) -^i / 

When {Aj} = {ni}, Eq. (48) gives the Gaudin formula for the norm of a Bethe state 

[21, 52] 

MM M . , _ X I ^ 

(Oin^(^.)n^(^.)|0) = «i^l^''^ n ''Iwa ^\7^ detMg-({A}). (50) 
jJi jJi abti smh(Aa-A6) 



where H' is 

H'M\}) = -S,t"' " ' 



a=l 



i^(A,-Afe), (51) 



and 



= ■ "^y't. ^ • (52) 
^ ' smh(A + ?7) smh(A - r?) ^ ' 

Notice that for the following manipulations, it is fundamental that the set of numbers 
{/ij} could not be solution of some Bethe equations. 

3.2. Reduced density matrix 

Let us consider a given Bethe state | {Aj} ), and let us select a block of length f as a 
subsystem of the spin chain. Every element of the reduced density matrix of these £ 
contiguous spins can be written as 

where the indices e can have the values {+, — } and the matrices E"^'^ are 

Once we know the reduced density matrix, any multi-point correlation function built 
within the £ spins can be found by considering the appropriate linear combinations. 
The most general object we need is 
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that has been obtained in Ref. [23] 



where 

Before reporting the main result of [23], we have to define the following two sets of 
indices 

a+ = {j:l<j< £, e, = +} , (56) 
a- = {j:l<j< £, e'j = -} . (57) 

We denote with rf+ (rf~) the dimension of the set a+ («")• For each j e it is 
necessary to define a set aj (if j G a~) and a set a'^ (if j € a"*") such that 

1 < < M + j, Uj eAj, l<a'j <M + j, a'j e A'^. 
where we introduced 

Aj= {b : 1 < b < M + e, b^ ak,a'k, k < j}, (58) 

A'j = {b:l<b<M + e, b^ak,k< j, b^Uk, k< j}. (59) 
Now we need only the redefinition 

{Xk}^{\k,^i...^i}, (60) 

to write [23] 

M 

( I n C(Afe) T,,,,; (Am+i) . . . (Am+£) = (61) 

k=l 

J2 G{,^,a^}(Ai,...,AM+^)(0| n C'(^b)' 

where 

M+i-i 
n sinh(A„^. -Xb + T]) 

b = l 



G{a,,a^}(Ai,...,AM+«) = Yi d{Xa,)- j^^^ 

jea- sinh(Aa^. - A6) 



b=l 

M+j-1 

n sinh(Afe - Xa'. + rf) 

6=1 ^ 

n 
• (62) 

jea+ n sinh(A6-A„0 

b=i ' 



An important simplification comes from the relation (42) which allows 
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Now it is important to know how many terms are involved in the summation in 
Eq. (61). By simple counting we get 



Ml 



-d- +a+ -i+1) , (63) 



(M-rf-)!ll 

that in the limit of large M behaves as . 

We stress now one of the main features of this approach for the calculation of the 
reduced density matrix. The computational resources we need for the algorithm grow 
exponentially with £, so it would be comparable to exact diagonalization, but at fixed 
i they only grows algebraically with M (but with a power equal to £). Thus we can 
expect that for relative small £ we can calculate the reduced density matrix for very 
large systems, while exact diagonalization can work with at most about 30 spins. 

Thanks to the invariance under permutations of the set {Ai,...,Afc} in the 
Slavnov formula, the number of determinants we need to calculate can be reduced 
to 

d+ f d+ 



f ''^ r ' 

h<h<---<3i '-fe=l 



M 

d~ +d+ -i 



(64) 



1=1 ^ jl<j2<---<ji 

We exploited this symmetry to reduce the computational effort. 
To proceed further, we define the set A 

A={{ai,...,aa-+d+)}, (65) 

whose elements have the property that a, ^ aj and 

OiG {!,..., M} if l<i<d-, (66) 

aiG{l,...,M + af} if d' < i < d' + d+ , (67) 

that allows to rewrite the summation in Eq. (61) as a sum over the elements of the 
set A 

E = E • (68) 

{aj,a'.} {ai,...,a^-^^+}eA 

Using the definition of Aj, we can write 

(0| n C{X,) = {K,,K,...,Xi^_^,ij,,...,^,J (69) 

66Ai + i 

where 1 < ii < M and 1 < ji < I. Wc stress that in general I. 

In order to calculate (55) we have to take the scalar product between (69) and 
the state 

ns(A,)|0). (70) 
Using again (44) we can rearrange the set of Aj 

{-^il ) '^j2 5 • • • I "^iu-n 1 ^ki ) • • • ) -^fe„ } ) (71) 

where 1 < ki < M. The important point is that the first M — n rapidities Aj are the 
same as in Eq. (69) and give a trivial contribution to the scalar product. 
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3.3. Homogeneous limit 

The tricky task in dealing with algebraic Bethe Ansatz is to take the homogeneous 
Umit in Eq. (61). In order to perform this limit, wc consider first the Slavnov formula. 
We remind that in general we have to do the scalar product between the two states 

I {ft} ) = \{Ki,K2^ • • • 7 ^ijv-„'?ii) • • • »Ci„}) ) 

I {^i} ) = \{^ii,K2J ■ • ■ ) ^iN-n^^ki, ■ ■ ■ 1 •^fe„}) 1 

where < n < i and ig, kg can take values in the interval [1, M] while jg € [1, We 
have to specialize Eq. (48) to this case. It is easy to show that 

M M 

sinh(r?)^ n n sinh(Ai - m,- + v) 



({AO Km,}) 



j=i i=i 



sinh(Ai — Xj) sinh(/Xj — /Ltj) 



M 

n n 

i=l j<'i 

where we introduced the matrix T whose elements are 



detT, 



j < M -n, 
j > M -n. 



(72) 



(73) 



sinh(Aa - fJ-b) sinh(Aa - fib + v) 
Recalling that for j > M — n we have /ij = , it follows that in the homogeneous 
limit we have a matrix whose last n rows are equal, so the determinant is zero. This 

M 

is compensated by the prefactor sinh(/ij — A*i) in Eq. (72) that is vanishing. 

i=l j<i 

To obtain the finite result of this limiting procedure let us define 

Mx) = /(A„ x) = sinh(A,-f-a;)sinh(Ai+3-x) ' ^^'^^ 

For an arbitrary value of n (ignoring for the moment the minus signs) we have to take 
the determinant of the matrix 

/ regular \ 

terms 



/i(ejj /2(eji) 



(75) 



where 



(76) 



(77) 



Let us now consider the Taylor expansion of the /(a, x) around x = 

f{a, x) = f{a, 0) + f{a, 0)x + ^f"{a, 0)x' + ... + ^/<"^(a, 0)x" + . . . . 

Gauss manipulations on the matrix above give the same result on each column except 
the index a which distinguishes the different columns. Therefore we can restrict to 
one column and construct the following matrix 

/ MO) rm^n ••• iT/fH0)e^V-- \ 



/i(o) /;(o). 



32 



Entanglement entropy of excited states 



25 



Since each column of the last matrix has the term fl'^\o) we can neglect it (we will 
restore it at the end of the manipulations) and consider the matrix 

1 e ^e? 



P 31 
I J 



V 1 



(78) 



By mean of rows manipulations it is possible to put the last matrix in a triangular 
form 



/ 1 





V 



5i(eji) 

52(eii,ei2) 






5!52(e|i,e^J 

1 _ (Jl J2 Jl \ 

2l53(eji,ej2'S-3'' 



i9i(4) 

1 r, I J J \ 



\ 



(79) 



where we have for instance gi{x) = x, g2{x,y) = x — y and more complicated 
expressions for the other functions. In order to obtain the homogeneous limit in 
the Slavnov formula, we need one more step; since we know that the limit exists, we 
have to choose in a convenient way the variables Cj^ . One possible choice is 



f + eexp(^j) j>l 

that corresponds to = eexp(^^jfc). It is useful to consider the simpler case in 
which Cj. = Cj. Substituting in (78) we obtain that the matrix (79) has a simple form. 
Iiidccd it is easy to see that (79) becomes proportional to the identity matrix (up to 
the n-th order) Kl , with 

[exp(l^jfc)- 



K = det 



(80) 



In conclusion this means that to have the lowest order in e for the Slavnov determinant 
we can write the matrix (75) as 



Hi, 



K.f, 



U-M+n-l) 



j < M -n, 
j = M -n+1, 
(0) j > M -n+1. 



Moreover we have to consider the contribution given by 



n sinh(a-C,) = (-e) 



n(n-l)/2 



j>k 



n 

fc=i 

j>k 



(81) 



(82) 



This concludes the calculation of the homogeneous limit in the Slavnov formula. 

We can now ask how many terms is it possible to obtain with the algorithm 
developed so far. The answer can be given examining the function G in Eq. (61). If the 
function G has no poles, then the procedure just outlined works with no modification. 
Unfortunately this happens only in very few cases, for example for the first element 
of the reduced density matrix, that is the so called emptiness formation probability 



{^9\ n i(i-cTi)iVs) 



(83) 
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e 


Here 


DMRG 


1 


0.49999999999999 


0.5 


2 


0.17659666969479 


0.17659666969468 


3 


0.04110985506014 


0.04110985506012 


4 


0.00595577151455 


0.00595577151455 


5 


0.00050690054232 


0.00050690054232 


6 


0.00002367077112 


0.00002367077112 


7 


0.00000055351689 


0.00000055351689 



Table 1. Emptiness formation probability of a chain of 20 spins in the ground 
state for A = 0.5 (on the left). We compare our result (left) with those 
(numerically exact) from DMRG in Ref. [55]. 



(On passing it is worth mentioning that this clement can be computed in the 
thermodynamic limit [53, 54, 35], basically because of this simplification.) In this 
case Eq. (61) simplifies to 

M e M 

{Q\ n C{Xa) n D{^j) U B{Xa)\0} 

= MW) M M ' (84) 

{0\UC{Xa)UB{Xa)\0) 
a=l a=l 

that can be written as 

M £ M+1 M+2 

( 1 n ^(^'^) n ^(^^+^) = E E • • • (85) 

k=l j=l ai = l 02 = 1 

027^ai 

M+i M+i 

... Yl Ga,...uMl--->^M+i){0\ n ^(^'^)' 

ai=l k = l 

"iT^oi oi-l fe^oi,...,o; 



where G is 



n sinh(Aa^. - Xb + r]) 



G„,...„,(Ai,... /aM+£) = n^(^«.)^^T^7 • (86) 



M+i 

^■=^ n sinh(A„, - Xb) 



b=l 



By definition G cannot diverge, then all the machinery developed so far is enough 
to compute t{£). In Table 3.3 we report some results for the emptiness formation 
probability for a chain of length L = 20 at A = 0.5 in the ground state. The agreement 
with DMRG data is perfect (taking into account the numerical rounding-off). 

There is another class of elements of pi accessible without further manipulations. 
In Eq. (61) the term that is easily manipulated is 

M+j 

H sinh(Ab - A„p , (87) 

b = l 

thus the only class with no poles is when in the sets Aj aj . . . we have 

a+ = {1} , (88) 
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corresponding to the elements -Pi o o") which is one particular column of the reduced 
density matrix. For the other 2^ — 1 columns we need still further manipulations. 

3.4- A last trick for the general case 

The problems in the general case arise from the divergencies of the term (87). Let us 
start with some preliminary observations. First, it is important to know the maximum 
degree of the poles in (87). Given a term of the summation in (61), the order of the 
pole is the order of the zero in 

M+3 

n J] sinh(A6-A„p, (89) 



that is given by 

E K^-^)' (90) 

where ffa'^ is the number of elements j € a'^ such that a'j > M. It is easy to max;imize 
the last expression to find 

d+ 

E'K+-i+i-^)' (91) 
1=1 

where the prime means that the summation is restricted to the i such that a'^+_^_^_^ — 
i > 0. 

We construct the general procedure to tackle the poles in (87). In the following 
we report our solution to the problem. However before proceeding we stress that this 
solution is somehow unsatisfactory, because we will end with and enormous sums of 
determinants for each element of pg. If we would have been able to find a "shortest" 
representation of the same elements, we could have been able to describe much larger 
£. Further developments in this direction would allow this method to be competitive 
even with DMRG [56] for the ground-state. 

Since the determinant in front of the pole is in general finite and the final result 
must be finite, all the coefficients multiplying each pole must sum to zero in (61). 
Furthermore this implies that we can ignore these terms (because we know in advance 
that they give zero) and concentrate on the important ones. To proceed, it is necessary 
to reshuffle the various terms in (61). Let us define 

M M 

n n sinh(Ai - iJ,j + rj) 

Q _ j = l 8=1 ^ .g^N 

n sinh(/ife - jjLj) n sinh(A/3 - A„) 

j>k a</3 

and 

f = detT. (93) 

We know that T ~ Jn the homogeneous limit. However here we have in 

general a pole of order n(n — l)/2 + q in G, then we have to expand both T and 
the nonsingular part of G up to the order n(n — l)/2 + q. For T, we developed the 
following procedure. Instead of doing the substitution (81), we put the higher orders 
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N 


Here 


Exact [49] 


Here 


Exact [49] 


27 


0.4130835714633 


0.4130835714633 


0.1879727171090 


0.1879727171090 


51 


0.4108297243638 


0.4108297243637 


0.1851632322689 


0.1851632322688 


101 


0.4101798729742 


0.4101798729745 


0.1843536264631 


0.1843536264633 


151 


0.4100571750358 


0.4100571750361 


0.1842007880727 


0.1842007880729 


201 


0.4100139161598 


0.4100139161591 


0.1841469044722 


0.1841469044717 



Table 2. Trp| (left) and Trp| (right) for the ground-state of A = 0.5 obtained 
here compared with the known exact results. 



in e up to n + q. Doing so we know that the determinant gives a polynomial in e with 
lowest degree is nid = n{n — l)/2, and we indicate with the highest degree. Thus 
the determinant is a polynomial of the form 



f = am,e"'' + -.. + aM,e^' (94) 

We can calculate all the coefficients ai numerically: it is enough to calculate the 
determinant in _D = Md — rrid different points and then to invert the linear system. 
Moreover if we choose the points in a smart way 

(27ri \ 
— fcj (95) 

the solution of the linear system is numerically trivial since the matrix of the system 
is unitary. Using this procedure we are able to calculate, in principle, all the elements 
of the reduced density matrix. In practice, our possibilities are limited by the size of 
the density matrix. Actually for small sizes we can go quite far and obtain the pf with 
three spins for chains with 200 spins, a task impossible with exact diagonalization. 

In Table 3.4 we show the quantity Trp^ for n = 2 and n = 3 for odd chains at 
A = 0.5 where we can compare with the exact results in Ref. [49]. The agreement 
is perfect and the small differences are due to the numerical rounding-off (we are 
summing order of 10^ elements in double precision 10~^^). 



3.5. Results: Entanglement entropy of excited states 

The main advantage of the method we have developed in the previous subsection 
is that we can exactly evaluate all elements of the reduced density matrix for any 
eigenstate of the XXZ chain. Compared to exact diagonalization we do not need to 
fully diagonalize the 2^ x 2^ matrix to find the eigenstates, we can just pick up our 
desired state by choosing the correct quantum numbers. As we have already pointed 
out, once the eigenstate has been chosen, the numerical complexity of the algorithm 
is only a power-law in TV (actually M, but for the most interesting states they are 
proportional), but the exponent grows linearly in I limiting the range of applicability 
of the method. If we would have been interested only in the ground-state properties, 
this method is less effective than DMRG or any method based on matrix product states 
[56]. In fact, these numerical methods require very little numerical effort to get the 
spectrum of the pe at machine precision for the system sizes that are accessible to us. 
However it is hard, if not impossible, to calculate the entanglement properties of highly 
excited states with DMRG. Thus our method, based on algebraic Bethe Ansatz, is by 
far the most effective available. We checked that our algorithm reproduces the known 
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Figure 8. Entanglement entropy of the excited states of the XXZ spin-chain 
for A = 10~^, 0.1, 0.3, 0.5 with N = 24 plotted against the logarithm of the 
conformal distance. Left: States that in the fermionic description for A = have 
two discontinuities. The slope agrees with effective central charge equal 1. Right: 
States that in the fermionic description have four discontinuities. The results are 
compatible with an effective central charge equal to 2. The bottom blue-line is 
the entropy of the ground-state at A = 0.3 shown for comparison. 



results for the ground-state for several different A, but we do not find instructive to 
report these results here. 

For the study of excited states, we consider spin-chains of length N = 24 
in the critical antiferromagnetic region (0 < A < 1) for four different values of 
A — 10^^,0.1,0.3,0.5. Using our algorithm we generate the full reduced density 
matrices with £ < 6 spins and from this we calculate the entanglement entropy (for 
the ground-state we know that already these small values of £ capture the asymptotic 
behavior [49, 55]). 



3.5.1. Bethe equations and classification of the states. The Bethe equations (43) 
can be re-casted in a form that is useful for numerical solutions and for a complete 
classification of the states. For practical reasons, we consider only the number of sites 
A'^ to be even. We recall that the Hilbert space separates in sector with defined number 
of reversed spins M (with respect to the reference state cf. Eq. (41)), that gives the 
total spin of the state in the z direction S'^^'^ — N/2 — M . Taking the logarithm of 
Eq. (43) and posing C, = arccos(A), we have 



at an 



tanh(Aj) 
tan(C/2) 



1 



M 

-V, 

fc=i 



at an 



tanh(Aj — A^) 
tan C 



N 



(96) 



Each set of distinct half-odd integer (integer) for M even (odd) numbers {It] (defined 
mod(A^)) specifies a set of rapidities, and therefore an eigenstate. For example, in the 
ground state these numbers take the values 



r(0) 



Af + 1 



+ i , j 



(97) 



This ground state can be interpreted as the spinon vacuum. Spinous are the 
elementary excitations of the model. They have spin 1/2 and obey semionic exclusion 
statistics (see e.g. [57] for a simple introduction to these excitations). Excited states 
have a defined number of up-spinon and down-spinon n_ . The total number of 
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spinons n+ + n_ < A'' is even when N is even, while n+ — n_ = N — 2M = 25j'^^. 
(Actually for the interacting model different values of S^qj, are possible when the state 
has some higher strings, which are non-dispersive. This discussion is too technical 
for the goals of this manuscript and we remand the interested reader to Ref. [43]). 
However, we will see that the spinon content is not the most important quantity for 
the entanglement entropy of excited states. 

Employing the property that the quantum numbers Ij are defined mod(iV), we 
can choose the allowed ones in the sets 



Only a subset of these numbers, bounded by a calculable Imax function of A and M 
(see again Ref. [43] for the technical details) provides real solutions for the rapidities 
A, and we limited our attention to these states. Fixed the parity of M (i.e. of S^qj,, 
since N is even), any state is defined by taking M numbers among the allowed ones 
in j(°dd) Qj. j(even)_ r^^iQ spinon content of the state then follows (see again [57]). 
Instead of using this standard Ij notation to indicate the states, following Ref. [58], 
we adopt a more complicated one that is useful to recover the fermionic description 
of the XX model when A = (because we want to compare with the results in the 
previous section). We denote with the spinons with polarization up (down) 

and with O^^"-* the empty positions that can be occupied by spinons with up (down) 
polarization. We indicate with an exponent the number of consecutive symbols, for 
instance •^+0^^"'' stands for 



To each sequence of N of these symbols • and G we can associate a single 

configuration of 7^. For example, the sequence 0^^^#'^^0^^^ has two up-spinons 
and no down ones; it follows that we need M = N/2 — (n+ — n^) = 11 quantum 
numbers; the state is fixed by taking from the set 1°'^'^ the last 11 numbers (see again 
Ref. [58] for more details). The advantage of this maybe not really intuitive notation 
is that the rule to recover the fermionic description is very easy: given the sequence 
one has to associate a fermion for every #~ or O'*' [58]. 

Once we have the set li for each state we are interested in, we use the Newton 
method to solve the Bethe equations. We limit ourself to the two-spinon and four- 
spinon sector of the spectrum (but we could easily consider other states). We also 
select states with real rapidities, to avoid problems with strings contributions, that 
however can be handled following Ref. [43]. 

3.5.2. Results. The main feature we want to check here is if the conformal scaling 
(2) with an effective central charge a is still valid for given excited states when we add 
the interaction A to the XX chain considered in the previous section. The prediction 
for the XX is based on the discontinuities of rfi{(t)) (cf. Eq. (19)). In order to predict 
the result at A ^ we exploit the mapping between the fermionic description and the 
spinonic one at A = 0. Once we have the fermionic picture associated to the state, 
we have to((/)) for A = 0. In order to check if the logarithmic scaling is obeyed, we 
plot Se, against Sf^ in Eq. (33), so that if the dependence is linear, the slope gives 
automatically the central charge a of the effective Hamiltonian. In Fig. 8 (left) we 




(98) 



•+•+0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 + 



(99) 
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Figure 9. Entanglement entropy for the two- and four-spinon states with 
A = 10~^, 0.1, 0.3, 0.5. Left; Summary of all the states we considered (for space 
problems, the legend shows only states at A = 0.3). Right: Independence of the 
leading term on the spinon polarization. We considered S^qj, = 0, 1, 2. The slope 
does not depend on the polarization. The bottom-red line is the ground-state at 
A = 0.3. 



display some states in the two-spinon and four-spinon sectors. We choose these states 
in such a way that in the Hmit A ^ 0, the corresponding fermionic structure has two 
discontinuities. For example, the state 0®~*^^0^"'^'^0^~ corresponds to the fermion 
representation | i^t^^i^)i having two discontinuities in m(ip). For A = 0, we know 
from the previous section, that all these states are described by Eq. (2) with effective 
central charge a = 1, as in the ground-state. Fig. 8 (left) provides a clear evidence 
that the asymptotic behavior of the entropy for i ^ 1 does not depend on A, at least 
in the considered range A g [0,0.5]. In the figure we also report the ground-state 
value for A = 0.3 for comparison. 

In Fig. 8 (right) we report the entropy for some states whose fermionic description 
contains four discontinuities. Again we can observe that the data support the 
logarithmic behavior in Eq. (2). The slope is different from the ground-state one, 
and indeed a naive fit (i.e. ignoring further corrections to the scaling that at ^ < 6 
are important) of the constant a gives a ^ 2.3, which is in agreement with the XX 
prediction a — 2. Moreover, the state 0^~0^'''#^^0^~''0^~ shows that the additive 
constant c'l in Eq. (2) depends dramatically on the details of the state (as we already 
know in the XX model). In Fig. 9 (right) we show the dependence on the spinon 
contribution of the additive constant. In Fig. 9 (left) we report the von Neumann 
entropy for all states and values of A we calculated. The changing in behavior for 
different numbers of discontinuities is clearly visible. In this figure we also report 
two (almost indistinguishable) states that have six discontinuities in the fermionic 
description and so are expected to have a = 3. There are strong crossover effects 
preventing us to extract clearly the value of a for such small subsystems, but the 
data are clearly in the right direction. This crossover is expected from the results 
for the XX model: when having 6 discontinuities in a chain of 24 spins, we expect 
approximately linear behavior in £ up to £* N/6 = 4, and in fact in the figure the 
crossover takes place around £ ^ 4. A quantitative understanding of this crossover 
(even for more excited states) requires larger values of £ and N that are not currently 
accessible to us. 

To conclude this section, we also report in Fig. 10 the data for log(Trp^) plotted 
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^ I \ I I I \ I \ u 

0.2 0.25 0.3 0.35 0.4 

l/41og(N/7i sin(7t£/N)) 

Figure 10. log(Trp|) against the logarithm of conformal distance. In the legend 
we only give the states for A = 0.3. 



against the logarithm of the conformal distance to check the conformal prediction [4] 

1 + n fN n£\ , 
-logTrp^^-^clog -sin- +c;, 100 
on \Tr 1\ J 

for n — 2. In fact, if the slope of all the previous curves can be interpreted as the 
central charge of some effective critical Hamiltonian having this state as a ground- 
state, not only the entanglement entropy should follow the conformal prediction (2), 
but also all Renyi entropies should scale according to Eq. (100). And in fact, as for 
Sf, the curves arrange in sectors with approximately similar slopes. Strong even-odd 
oscillations of the Renyi entropies prevent us from any reliable quantitative analysis, 
as it is the case in the ground-state [49] . Again it is visible the same structure observed 
for the von Neumann entropy. However, naive fits give reasonable estimations of the 
effective central charges for the two lowest sets, but the oscillations (combined with 
the crossover previously mentioned) spoil the result for the last set for which a = 3 is 
expect. 

The knowledge of the full reduced density matrix can be also used to calculate the 
entanglement spectrum (i.e. the distribution of its eigenvalues). However, because of 
the relative small values oil we can access, this is not enough to check recent conformal 
predictions for the spectrum [59]. 

4. Summary and discussions 

In this paper we considered the entanglement entropy of excited states in spin chains. 
We provided a full analytical study of the XY model in a transverse magnetic field. 
We found that the entanglement properties of the excited states depend strongly 
on the distribution of excitations above the ground state. To characterize them in 
the thermodynamic limit, we introduced the regularized characteristic function of 
excitation m(ip) e [—1,1]. The analytic properties of this function (or slight variation 
of it for critical systems) completely determine the entanglement in the scaling limit 
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N ^ £ ^ 1. When m{ip) ^ ±1 in a set of non-vanishing measure, we have that 
the entanglement entropy is extensive in the subset length (i.e. proportional to £). 
The analytic expression for such states is given by Eq. (17) as we proved by using 
the Szego lemma for block Toeplitz matrices. Oppositely when m{(f) = ±1 almost 
everywhere (as for the ground state) the entanglement entropy always follows the 
conformal scaling (2), even for non-critical systems (except in the ground-state when 
the area law holds). The pre- factor (that we indicate as a) is the central charge 
of a critical, local, translational invariant Hamiltonian, that we built explicitly. In 
the case of the XX model we proved this result rigorously via the Fisher-Hartwig 
conjecture. These logarithmic states have a finite-size scaling that is by construction 
the conformal one in Eq. (2). Oppositely the extensive states have very peculiar finite 
size scaling with slopes that changes according to the analytic properties of m{ip). We 
have been able to connect these features to the (non-)locality properties of an effective 
Hamiltonian that can be made local on a wrapped chain. 

We also considered the XXZ spin chain, that is solvable by Bethe Ansatz. We used 
the algebraic construction to calculate exactly the reduced density matrix for finite 
chains with I < Q. The method we developed is ideal to obtain the entanglement 
entropy of excited states. In fact, while numerical methods based on MPS like DMRG 
[56] are very effective for the ground-state, they usually work bad for highly excited 
ones. Our method instead treats on the same foot any eigenstate, that is specified by 
the quantum numbers related to the spinonic content of the state. This method has 
the numerical advantage that its complexity increases only in a polynomial way with 
N (while exact diagonalization is exponential). The drawback is that the complexity 
increases exponentially with £ and limited our study to £ < 6. We do not not know 
whether this is an intrinsic limit of the method, or if our representation of the reduced 
density matrix can be still drastically optimized to make the procedure more effective. 
The trickiest point in our derivation was to obtain the homogenous limit from the 
results in Ref. [22, 23]. If we would have been able to find a more effective way to 
perform this limit, the method we propose could have been as effective as DMRG. 
However, even if we could study only subsystem with £ < 6, we have been able to 
conclude that the main results obtained analytically for the XX model (at A = 0) 
remain valid when interaction is turned on. We showed in fact (making the proper 
mapping between spinonic and fermionic excitation at A = 0) that all the states 
that are logarithmic for A = maintain this property with the same prefactor and 
with a non-universal additive constant that depends very smoothly on A (as for the 
ground-state [55]). 

After this study, the characterization of the asymptotic block entanglement of 

excited states in these two chains is at an advanced level. Few unsolved problems are 
still present, especially for the XXZ chain, as e.g. the understanding of the string- 
states and the quantitative description of the crossover between linear and logarithmic 
behavior. However, the main question that still remains open is how general are these 
results. The fact that the division among extensive and logarithmic states is conserved 
when the interaction A is introduced, strongly suggests that this phenomenon should 
be expected for any local spin-chain, with a prefactor that can be predicted after 
that the relevant excitations have been identified. In fact, in the interacting system 
(especially for A not small) the excited states are complicated linear combinations of 
the free-particle ones, several degenerations are also removed by A. and it is unlikely 
that such result is only a coincidence. However, we do not have a general proof for 
this statement. It might be that for low-lying excited states the generalization of the 
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conformal methods of Refs. [4, 5] can give the logarithmic behavior of these states. 
Anyway, we have shown here that this property is not Hmited to low-lying states 
and so a more general proof would be desirable. Another very interesting question 
would be to understand how these results are affected by quenched disorder. For 
the ground-state it is known that only the prefactor of the logarithm is changed [60] , 
but the excitations in these systems are so different that major qualitative changes 
can take place. Finally, it should be possible to generalize the methods employed for 
ground-states of free lattice models in higher dimensional systems (as those reviewed 
in Ref. [2]) to the excited states of the same models. Some results in this direction, 
relevant for the physics of black holes are reported in Ref. [9] . 
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